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Summary Integral representations are derived for the parabolic 
cylinder functions U(a, x), V(a, x) and W(a, x) and their derivatives. 
The new integrals will be used in numerical algorithms based on 
quadrature. They follow from contour integrals in the complex plane, 
by using methods from asymptotic analysis (saddle point and steep- 
est descent methods), and are stable starting points for evaluating 
the functions U(a,x), V(a, x) and W(a,x) and their derivatives by 
quadrature rules. In particular, the new representations can be used 
for large parameter cases. Relations of the integral representations 
with uniform asymptotic expansions are also given. The algorithms 
will be given in a future paper. 
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1 Introduction 

The solutions of the differential equation 



are called parabolic cylinder functions and are entire functions of z. 
As in [1], Chapter 19, [7], and [9] we denote two standard solutions of 
(1.1) by U(a, z), V(a, z). Another notation found in the literature is 
D u (z) = U(—v — 2,z). Special cases are Hermite polynomials, error 
functions and Fresnel integrals. 




(1.1) 
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Values at the origin are given by 



U(a,0) = ^f ,U'(a,0)- ^ 

1 „4- 1 1 



2i a +ir{ 3 + l a) > V >-/ 2 i-| r( l + l a) ; 



(1.2) 



Then we have 

I7(a,z) = [/(a,0)yi(a,z) + [/'(a,0)y 2 (a,z), (1.3) 
F(a,z) = y(a,0)yi(a,z) + y'(a,0)y 2 (a,z), (1.4) 



where 



1,2 



(1.5) 



yi(a, z) = e 32 1F1 (-£a -^z 2 ) 
= e" l22 1 F 1 (^ + ii;^ 2 ), 

l z 2 , v 

y 2 (a,z)=ze 4 1F1 (-£a + |, §; -^z 2 J 
= ze~ 32 iF^^a + if;^ 2 ), 
and the confluent hyper geometric function is defined by 

n=0 V C M n - 

with (a) n = r(a + n)/r(a),n = 0, 1, 2, . . .. 

The functions yi(a, z) and y 2 (a, z) are the simplest even and odd 
solutions of (1.1) and the Wronskian of this pair is given by 

W[ yi (z),y 2 (z)} = yi(z)y' 2 (z) - y' 1 (z)y 2 {z) = 1. (1.7) 

/,From a numerical point of view, the pair {yi, y 2 } is not a satisfactory 
pair [4], because they have almost the same asymptotic behaviour at 
infinity. 

The behaviour of U (a, z) and V(a, z) is, for large positive z and 
z » \a\: 

U(a,z)=e-* z2 z- a -h\+0(z-i)], (Lg) 

V(a, z) = z^ [1 + O (z- 2 )} . 

Clearly, numerical computations of U(a, z) that are based on the 
representations in (1.3) and (1.4) should be done with great care, 
because of the loss of accuracy if z becomes large. Also, for large a 
these representations become useless. 
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The Wronskian relation between U(a,z) and V(a,z) reads: 



W[U(a,z),V(a,z)] = J2/ir. (1.9) 



W[U(a,z),U(a,-z)]=-f- T - (1.10) 

which shows that U (a, z) and F(a, z) are independent solutions of 
(1.1) for all values of a. Other relations are 

U (a, z) = 2 — — p- [V(a, -z)- sin ira V(a, z)\ , 

cos ira 1 (a + 77 ) 

r (o + i) 

V(a, z) = ^ [sin7ra U(a, z) + £7(a, —z)] . 

Equation (1.1) has two turning points at ±2yf^a. For real pa- 
rameters they become important if a is negative, and the asymptotic 
behaviour of the solutions of (1.1) as a — > —00 changes significantly 
if z crosses the turning points. At these points Airy functions are 
needed for describing the asymptotic behaviour. 

The purpose of this paper is to give integral representations of 
U (a, x) and V(a, x) for real values of a and x. We use integral rep- 
resentations from the literature and modify these by saddle point 
methods. In this way we obtain integrands that are non-oscillating, 
also for the case a < 0. In particular, we can use the new representa- 
tions for large parameter cases. In earlier papers [8] and [2] we have 
used these methods for obtaining stable integral representations for 
modified Bessel functions with pure imaginary order and for inhomo- 
geneous Airy functions (Scorer functions). 

We give relations of the integral representations with uniform 
asymptotic expansions, which are taken from [6] and [10]. We only 
give the expansions in terms of elementary functions. Uniform expan- 
sions in terms of Airy functions can be found in [6] , and a modified 
form in [10]. 

We also consider solutions W(a, ±x) of the differential equation 

W" + {^x 2 - aj W = 0, (1.12) 

a modified form of (1.1), again for real a and x. Properties of W(a, x) 
are given in §4, which can be found in [1] and [5]. 

In a future paper we give algorithms based on quadrature rules 
for evaluating the integral representations of U(a,x), V(a,x) and 
W(a,x). 
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In [10] numerical and asymptotic aspects of the parabolic cylinder 
functions have been discussed, and we refer to this paper frequently. 
The notation of certain quantities is also as in [10]. The asymptotic 
methods referred to in this paper (saddle point methods) can be found 
in [7] and [11]. For an overview of the numerical aspects and software 
for the parabolic cylinder functions we refer to [3]. 



2 Integral representations for a > 

We derive integral representations for U(a, x) and U(a, —x). The com- 
putation of V(a, x) for a > can be based on the second relation in 
(1.11). For a > the functions U(a,x) and U(a,—x) have a non- 
vanishing Wronskian relation (see (1.10)), and moreover, these func- 
tions constitute a numerically satisfactory pair of solutions of (1.1). 



2.1 The case x>0 

We take the integral (see[l], formula 19.5.4) 

U{a,x) = ^= I e-**h*s-*L (2.1) 

where C is a vertical line on which 3fts > 0. On C we have —^tt< 
phs < tj7t, and the many- valued function s~ a_1 / 2 assumes its princi- 
pal value. The transformations 

x = 2t^/a, s = y/aw (2.2) 



give 



where 



!,(.,«) = £*^g± (2.3) 
z\/2vr Jc V w 



(f){w) = ± w 2 -2tw-\nw. (2.4) 

The saddle points follow from solving 

,// \ w 2 — 2tw — 1 , . 

4>'{w) = = 0, (2.5) 

giving saddle points at t ± Vt 2 + 1- We take for the path C in (2.3) 
the vertical through the positive saddle point 



w 



= t + Vt 2 + l. (2.6) 
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At this saddle point C coincides with the steepest descent path trough 
wo. The complete steepest descent path follows from solving ^[^(iu)] = 
Q[4>(wo)]. In the present case 9?[(/>(w>o)] = and we obtain for the sad- 
dle point contour the equation 

±r 2 sin26» - 2trsm6 - 6 = 0, where w = re ie , (2.7) 



which can be solved for r = r{6): 

t + ^t 2 + 6 cot e 

v — 

cos 9 



-§7T < 6 < ±7T. 



(2.8) 




Figure 1. Steepest descent contour for the integral in (2.3). 
Then (2.3) can be written as 

l x 2 +a ^,( W0 ) i_i 1 

U(a, x) = -j= J ^ e"Mg{0) M, (2.9) 

where 

^(6») = $l[(j)(w)-(j)(wo)] = ±r 2 cos26-2trcos6-lnr-cf)(wo), (2.10) 



and 



9(0) = 9 



1 dw 



^ (2cosfl+ l)r 2 -2tr + l 
VrcosiflV^ 2 + 0co7^ 



(2.11) 
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The quantity £ denned by 

£= |[t V / t 2 + l + ln(t+ Vt 2 + 1)], 



(2.12) 



is used in the asymptotic representation of U(a,x) in this case; see 
[10], formula (2.29). We have 

\x 2 + a(f>{w ) = a[±-Wt 2 + l-\n(t + Vt 2 + 1)] =a(i-2|). (2.13) 



This gives 



[/ (a, x) 



2< /.Itt 



a* e 



2tt7( 



- / e a ^g(e)de, (2.14) 



where 



7 (a) = e~2 a a2 a . (2.15) 
For the derivative U'(a,x) we can start from (2.1), and we have 



U'(a,x) = 
This can be written as 



I-r 2 



W2tt Jc 2 



(2.16) 



U'(a,x) = - 



a* e 



-2a£ 



27T7(a) 



e a ' m h{9)d9, (2.17) 



where 



1 du> 



r 3 - tr 2 (2 cos - 1) + r(2i 2 + 1 + 2 cos 9) - t 
4^¥cosUVt 2 + 9cot9~ ' 



(2.18) 



2.2 The case x < 

This case can be done by using the representation of the previous 
section. However, when t is a large negative number, the saddle point 
wq defined in (2.6) is close to origin, at which point the integrand 
of (2.3) is singular. As a consequence, the functions ip(9) and g{6) 
in (2.9) have singularities close to the origin 9 = when t is a large 
negative number. 
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In the present case we start with the well-known integral repre- 
sentation (see [1], formula 19.5.3) 

-lx 2 

e f 00 a- A -ls 2 +xs , 

U(a,-x) = — 1- / s 5 e 2 ds, a > -±. (2.19) 

There are no oscillations, but it is convenient to transform the integral 
in such a way that the saddle point is at the origin and a suitable 
normalization is obtained. The transformations (2.2) give 



(2.20) 

where 4>(w) is given in (2.4). The positive saddle point wq is as in 
(2.6). We transform this point to the origin by writing w = wq(1 + u), 
which gives 

u{a , , x) _ '1 r ,-*<.> * , (M1) 

I (a + i) J-l y/l + u 

where we have used (2.13), 7(a) is defined in (2.15), and 

ijj(u) = \wlu 2 + u -ln(l + u). (2.22) 

For the derivative we have 

r(a + 5 ) (2.23) 



du 



u 



To avoid numerical cancellation for small values of u in the com- 
putation of ip(u) defined in (2.22), a specific code is needed for the 
evaluation of ln(l + u) — u. 



2.3 A Wronskian for the integrals 

When checking the numerical algorithms the Wronskian relations in 
(1.9) and (1.10) can be used. When the parameters are large it is more 
convenient to use a Wronskian relation that is based on the integrals 
derived in the section. This gives a better control of the errors that 
occur in the quadrature rules, because large and small factors are not 
present in the integrals. 
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We write (see (2.14), (2.17), (2.21), and (2.23), respectively) 

\ p -2a£ 

U(a,x)= I(a,x), (2.24) 

V27T7(a) 

at e - 2a « 

U'(a, x) = - I d (a, x), (2.25) 
V27r7(a) 

f/(a, -x) = V ° ^ / J(a, x), (2.26) 

J (a + 2 J 

U(a,-x) = v Jd(a,x). (2.27) 

A° + i) 

Then the relation for the integrals reads 

27T 

I(a,x) J d (a,x) + I d (a,x) J(a,x) = — — . (2.28) 



2.4 Uniform asymptotic expansions for a > 

The quantities related with the integrals are closely related with the 
uniform asymptotic expansions given in formulas (2.29), (2.33) and 
(2.34) of [10]. We have 

where F^(t), G^(t), Pfi(t), and Q^t) are supplied with asymptotic 
expansions that have a double asymptotic property: one of the pa- 
rameters o or t (or both) should be large. Recurrence relations for 
the coefficients of the expansions are given in [10]. 



3 Integral representations for a < 

We give integral representations for U(—a,x) and V(—a,x), with 
a > 0, and we consider three x— intervals. Let t = x/(2y/a). The dif- 
ferential equation (1.1) becomes for U(—a,2t^/a) and V{— a, 2t v / a) 
in terms of t 

^|-4a 2 (t 2 -l) y = 0, (3.1) 

which has turning points at t = ±1. Consequently, we consider the 
intervals t < —1, \t\ < 1 and t > 1. We start with the middle interval, 
where the oscillations occur. 
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3.1 The case -1 < t < 1 

We consider the integral 

r°° i9-i 

Y (a,x)= e~2 s +xls s a -2ds, > -\. (3.2) 
Jo 

Using (2.19), we see that 

Y{a, x) = r{a + l)e~^ x2 U(a, -ix). (3.3) 

We also have 

yi e - ¥ ia + \™U(a, -ix) = U{-a,x)/r{a+\) + iV{-a,x). (3.4) 

This follows from using the initial values in (1.2) and those of Y(a, x). 
It also follows from the relations in (1.11) and 19.4.6 in [1]. 
Hence, 

y| e -3^ a +3^ e ^V(o, x) = U(-a, x) + ir(a + \)V(-a, x). (3.5) 

We see that the single integral (3.2) produces U(— a, x) and V(—a,x) 
by taking real and imaginary parts. 

We proceed with Y(a,x), and the transformations x = 2 v / at, s = 
yfaw give 

Y(a,x) = a a / 2+1 / 4 r e -^^=, (3.6) 
Jo V w 

where 

(f>(w) = \w 2 - 2itw - Inw. (3.7) 
We consider a path through the saddle point 



w + = it + y/l - t 2 . (3.8) 

We have 

</>(w+) = \ + t 2 + 2i (ri - \tt^ , f]=\ (arccost - ty/l - t 2 ) , (3.9) 

where arccost has values in [0,7r] for t € [—1,1]- 

The path of steepest descent starts at tt; = 0, runs through 
and terminates at +oo; see Figure 2. The path follows from solving 
the equation 

3?0H = $Kf>(w+), (3.10) 
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that is, from solving 

\r 2 sin 20 - 2tr cos 9 - 9 - 2r] + \ir = 0, (3.11) 

where w = re ld . The solution of (3.11) reads 

t cos 9 + a^jt 2 cos 2 9 + sin 9 cos 9(9 + 2rj - \tt) 

sin 9 cos ' (3.12) 

< 9 < O , 

where 0o = —2?] + ^tt; the square root is non-negative. The number 
a equals —1 when phw; + < < 0q, and +1 when < < phw + . 
Observe that phw; + = ^7r — arccost. When = 0o, we have r = 0; 
when = phu> + , we have r = 1, and = gives r = oo. For £ = the 
path coincides with the positive real axis. When t < the paths are 
in the lower half plane, and follow from those for t > by symmetry. 



1.5 








"~ ? = 0.9 


1.0 




^f_____ 




? = 0.5 


0.5 




w + 

. — • 




; = 0.1 




' i 

2.5 



Figure 2. Steepest descent contours for the integral in (3.6) 
for t = 0.1,0.5,0.9. 

A simple approximation of the path is given by (we write w 

u + iv) 

ut(l + u + ) 



u + u 2 . 



u+ = VI - t 2 . (3.13) 



This path runs through the point w + = u + + it, and has the same 
slope at this point as the exact steepest descent path, that is, dv/du = 
tj (1 + u+) for u = «+. 

For t = 1 the steepest descent path runs from the origin to w+ = i 
along the imaginary axis, and from i to 2i + oo. For more details on 
the case t > 1 we refer to § 3.2. 
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Integrating (3.6) with respect to 6 we obtain 
U(-a, x)+ ir(a + \)V{-a, x) = 

= ^1 a3 7 (a)e- i ( 2a?? -2 7r ) e^i/(«) <W, (3 ' 14) 

where 7(a) is defined in (2.15), 

^(0) = ir 2 cos2fl + 2trsinfl-lnr-t 2 - \, (3.15) 

and 

p ie / 2 I dr \ 

g{6) = _(_+ir). (3.16) 



We write the representations for U (—a, x) and V(— a, x) in real 
form, with trigonometric functions that correspond with those in [6] 
and [10]. We first write 

g(8)= gi (9)-ig 2 (e), (3.17) 
where gj(9),j = 1,2, are real. That is, by (3.16), 

do yjr 

cos(±#) drsin(±0) 

1BW = r _^i + __^i. (3.19) 

Then we have 

[/(-a, x) = y|a3 7 (a) [(sin A) Gi + (cos A) G 2 ] , (3.20) 

V(-a,x) = [Ij—^j- [(cosA)Gi - (sinA)G 2 ] , (3.21) 
V J l a ~r 

where 

A = 2ar/ + ivr (3.22) 

and 

G 3 = f e °e- a ^ gj (e)de, j = l,2. (3.23) 

JO 

For the derivatives we find, using (3.2) and (3.5), 
U'(-a, x) +%r(a + V'(-a, s) = 



2 a-^{a)e-^ 2ar i--^ / e~ a ^h{9)de. 
w Jo 



(3.24) 
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where 



h{6) = (t + iw)g{9) = h^O) - ih 2 {6). 
That is, by (3.25), (3.18) and (3.19), 



(3.25) 



h^O) = (t - r sin 6)g 1 (9) + r cos 9g 2 (0), 
h 2 {9) =rcos9 gi (9) - (t - r sin 9) g 2 {9). 



(3.26) 
(3.27) 



Then we have 



U'(-a, x) = J I a4 7 ( a ) [(sin A) H x + (cos A) H 2 ] , 



(3.28) 




(3.29) 



where A is given in Eq. (3.22) and 



Hj = r e- a ^hj(9) dd, j = 1, 2. 



Jo 



(3.30) 



3.1.1 A Wronskian for the integrals By using (1.9) and the integrals 
in (3.23) and (3.30), we find the relation 



r{a + \) = V2^ 1 2 {a)r*{a + \), r*(a + \) = 1 + 0(l/a), (3.32) 



as a — > oo. Hence, for large a, the right-hand side in (3.31) is of order 
^[1 + 0(l/a)] (see also formula (3.28) in [10]). The relation in (3.31) 
can be used for testing the numerical algorithms. 

3.1.2 Uniform asymptotic expansions —1 < t < 1 The relationship 
of the integrals Gj,Hj with uniform expansions follows from (2.23), 
(2.24) and (2.27) of [10]. These expansions are the same as in [6]. On 
the other hand, we can derive modified expansions (a main topic in 
[10]), by using (2.29) and (2.33) of that reference. From (3.4) and by 
changing t to —it in(2.29) of [10], we obtain 




(3.31) 



where 





(3.33) 
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U'(-a, x) +ir(a+ ^ V'(-a, x) = 

(3.34) 

-r* (a + ±) V2a3 7 (a)e-i^- 2ia "(l - t 2 )i G^-it), 
where F^—it) and G^(-it) have the asymptotic expansions 

^H*)~£(-ir^ ^Ht)~g ( -ir|p, (3.35) 

as a — > oo, uniformly for t & [—1 + 5,1 — 5]. The quantity r* is defined 
by 

T * = - 1 -(^JL= + l). (3.36) 



2 vyr^ 

The polynomials <j) s and ^ s are given in (2.11) and (2.16) of [10], with 
recursion relations. The first fraction at the right-hand sides of (3.33) 
and (3.34) has the asymptotic estimate 1 + 0(l/a) (see also formula 
(3.28) in [10]). 



3.2 The case t>l 

We use the the integral for Y(a, x) given in (3.6) with (p(w) given in 
(3.7). The saddle points are now purely imaginary: 



We have 



w-=it- iVt 2 - 1, w + = it + iVt 2 - 1. (3.37) 
4>{w±) = Jr + - Jrvri ± 2£, 



£ = \ [tVW^T - In (t + V^l) 



(3.38) 



The quantity £ is also is used in the asymptotic representation of 
U(—a,x) for this case; see [6] and [10]. 

The path of steepest descent starts at w = 0, runs through w- 
and w + on the positive imaginary axis, and from w + to +oo. The 
path from w + to +co follows from solving the equation 

Q(f>( w ) = -Ivr, (3.39) 

that is, from solving 



\r 2 sin 29 - 2tr cos 6 - 6 + ±tt = 0, 



(3.40) 
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where w = re ld . The solution of (3.40) reads 
t cos 6 + y^ 2 cos 2 6 + sin cos 6(6- \n) 



r = 



, < 6 < W. (3.41) 



sin cos 6 

The square root is positive, unless when t = 1 and = \-k 
We obtain 



C/(-a, s) + iF{a + l/2)V(-a, x) = a* 7 (a) 



e- a *Wg(6) d6 + ie 2 < f + e a ^ d ' 
Jo 



(3.42) 



where <p(v) = 4>{w^) — <j)(iv), g(6) as in (3.16), and ip(6) = 4>(w) — 
4>(w+) = dt[(j)(w) — <f>(w+)], w = re %e , now with r defined in (3.41). 
Explicitly, 

(j)(v) = \v 2 — 2tv + Inv — \r 2 _ + 2tr_ — lnr_, 
$(6) = \r 2 cos 26 + 2ir sin 6 - In r - \ - t 2 - 2£. 

where 

r_ = t - Vt 2 - 1. (3.44) 

Considering the real and imaginary parts on both sides of (3.42), 
we see that for V(—a, x) we need the v— integral with the dominant 
factor e 2a ^ and part of the 6— integral. When t ~ 1 (in fact, when a£ 
is small) both integrals are of the same asymptotic importance. The 
dominant saddle point in the v— integral is r_; in the 6— integral the 
dominant point is the upper limit. 

When we write 

e** i g(e)=g 1 (9)+ig 2 (9), (3.45) 
where g\{6) and g2(0) are real, we have 

U(-a, x) = yi a i 7 (a)e- 2 ^ G 1 (3.46) 

and 

V(-a, ,) = Mai^f^ (e-« C 2 + ft) . (3.47) 

where (for j = 1,2) 

Gj = ( h e-ofWgde) d6, G 3 = f + e a ^%. (3.48) 
Jo Jo V v 
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U'(-a,x) = ./jfa! 7 (a)e- 2< *#i, 



and 



V'(-a,x) = J^ 



7(a)e 



2a§ 



r 



(a+l) 



e~ 4 < F 2 + 



where (for j = 1,2) 



Hi 



[ h \-^hAe)de, h 3 = r + e *)((-4. 

jo ' jo v« 



(3.49) 



For the derivatives we have 
U'(-a, x) + ir(a + \)V'{-a, x) = fia^a) 

x [ h \- a ^h(e)d9 + ie 2 < ( r+ e ^)(t-v)^= 
Jo Jo Vv 

where h{9) = (t + iw)g{9). When we write 

e \™h{9) = hi(G) + ih 2 {9), (3.50) 

where h\{9) and h 2 {9) are real, we have 



(3.51) 



(3.52) 



(3.53) 



3.2.1 A Wronskian for the integrals By using the Wronskian relation 
in (1.9) and the integrals in (3.48) and (3.53), we obtain (cf. (3.31)) 



e-*<(Gifli - ff,G 2 ) + (Giffj - ffiG 3 ) 



a 7 2 (a) a 



(3.54) 



The relation in (3.54) can be used for testing the numerical algo- 
rithms. 



3.2.2 Uniform asymptotic expansions for t > 1 We give the rela- 
tionship of the integrals with the uniform expansions given in (2.9), 



Integral Representations for Parabolic Cylinder Functions 



17 



(2.14), (2.18) and (2.29) of [10]. We have 



Gi 



2^/E{t 2 -1)3 



Hi = - 



v^F(t 2 -l)4 
2v^ 



r{a + \)e a a~ a '^ 



(3.55) 



e- 4 < G 2 + G 3 



V2(t 2 - 1)3 



e- 4 < tf 2 + tf 3 



(t 2 -l)i 



v 7 ^ 



where F^t), G M (i), P^t), and Q M (i) are supplied with asymptotic 
expansions that have a double asymptotic property: one of the pa- 
rameters o or f (or both) should be large; t > 1 + 5. Recurrence 
relations for the coefficients of the expansions are given in [10]. 

3.3 The case t ~ 1 

For t ~ 1 the contours used in §3.1 becomes less suitable for nu- 
merical quadrature. For example, we see in Figure 2 that the saddle 
point w + approaches the imaginary unit when t | 1, and that the 
path becomes non-smooth when t = 1. For numerical calculations 
we may consider uniform Airy- type asymptotic expansions if t ~ 1, 
and we will investigate later if this is indeed the best approach. But 
we also investigate if a modified contour can be used for numerical 
quadrature. 

We use for t ~ 1 the representation of Y(a,x) in (3.6). We write 
w = u + iv, and integrate with respect to v along the line segment 
from the origin to w+, and then along the horizontal path from w + 
to W+ + oo with respect to u. In the first integral we substitute v = 
t(l — p), and integrate with respect to p. Observe that for t > 1 the 
point w+ is on the imaginary axis, and for this case no difficulties 
arise when t ~ 1, because the path is already split up into two parts; 
see § 3.2. 

It is not difficult to verify that the representations in (3.20), (3.21), 
(3.28) and (3.29) can be obtained, with Gj,Hj replaced with Gj,Hj 
(j = 1,2), where 
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Gj = /V^W gf\p)dp 

Jo 

+ / , T gf\u)du, 



1 e -a4 1} (j>) hf\p)d P 

roc p -aipr\u) , 

+ / ; ~hf (u) d«, 



where 

4 1] (P) = I^U - 2t 2 ) - P - ln(l - p), 

4 2 )(u) = iu + nVl-t 2 - | In (\ + 2uy/(l -t 2 + 



(3.56) 



(3.57) 



Y>1 j (p) =p 2 Wi-t 2 , 

(2) Ut 

w) (u) = arctan , — tu, 



9i\p) ~ = e* iT - a **i' ) V), t = arcsini, 

( 2 ) . ( 2 ) -^r-a#f ) W+^ arctan — U>.-.)8) 



h<?\ P ) = tpg?\p) + y/l=P(l - P)^ 1} (P), 
h { ?(p) = tpg^ip) - VT^(1 - p)g^(p), 
h?\u) = (u+VT^)g { 2 2) (u), 
hf\u) = -{u + ^l-t 2 )gf\u). 



(3.59) 



3.4 The case t < -1 

We can repeat the analysis, starting with (3.2) with x < 0, but do not 
need new integral representations, algorithms or uniform asymptotic 
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expansions for this case. For U(a, x) we can use the second relation 
in (1.11), and for V(a,x) the first relation. 

When the parameter a is large these relations have to be used 
with care, because gamma functions with large negative arguments 
occur. It is better to use the quantities Gj,Hj introduced in §3.2.1. 
In the computer code these quantities will be given as output from 
the case t > 1. 

We have 

U(—a,—x) = y^a47(a) je~ 2a ^ [cosvraG^ + sin7raGi] 

+e 2af costo G 3 | , 

U'(-a, -x) = -y^|a3 7(a) je~ 2a5 [cos7ra# 2 + sin7ra#i] 



+e 2af cosvra# 3 } , 



V(-a, -x) = ,/f ° 47( °i\ [ COS7raG, i " sinvraGs] 



2 



e s sm ira ( 



F'(-a,-x) = _ /I Q 4 7(a) r _2a€ [ COSTOi?1 _ s i n7ra tf 2 ] 

-e 2af sin to #3}. 



(3.60) 



4 The W- function 

In this section solutions of equation 

g + (ix" -„), = * (4.1) 

are considered, again for real a and x. For a < the solutions oscillate 
on the real x— axis; for a > there are turning points at ±2y/a, and 
the oscillations occur outside the interval [—2y/a,2y/a\. From quan- 
tum mechanics we know that (4.1) is the equation for propagation 
through a potential barrier. 



4-1 The standard solutions 

We consider solutions W(a,x) and W(a, — x); these form a numeri- 
cally satisfactory pair for —00 < x < 00; see [5]. The function W(a, x) 
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has the initial values (see [1], p. 692) 



W(a,0) = 2~4 



W'(a,0) = -Ti 



r{l + \ia) 



(4.2) 



The Wronskian of W(a, x) and W(a, —x) is 

W[W(a, x), W(a, -x)] = 1. (4.3) 

Power series expansions are 

W(a,a:) = W(a,0)w 1 (a,x) + W'(a,0)w 2 (a,x), (4.4) 

where ioi(a, x) and 102(0, x) are the even and odd solutions of (4.1). 
We have 



„2n 



wi{a,x) = ^a n {a)—— v 

n=0 ^ ' 



X 



2n+l 



w 2 (a,x) = £A.(o) (2w + 1) , 



, (4.5) 



where a n (a), f3 n (a) satisfy the recursion 

Q n+2 = aa n+ i - \{n + l)(2n + 1) a n , 
Pn+2 = ap n+1 - \{n + l)(2n + 3) /? n , 
Q (a) = 1, ai(a) = a, (3 (a) = 1, (5i{a) = a. 

The relation with the function U(a, x) reads 

-±-W(a, x) + z^/Fa) W(a, -x) 



(4.6) 



= v^ei^+^f/ (m, xe-^/ 4 ) , 



(4.7) 



which follows from using the initial values of the functions, but also 
from [1] [19.17.6 and 19.17.9]. The quantities k(a) and p(a) are given 
by 

k(a) = Vl + e 2 ™ -e na = , } , (4.8) 



and 



/9(a) = |vr + ±0 2 (a), <fc(a) = phT (± + ia) ; (4.9) 
the branch is defined by 02(0) = and by continuity elsewhere. 
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Because we assume that a and x, and hence W(a,±x), are real, 
we have, using (4.7), that 



W{a,x) = ^2k(a)e^ a ^ [e ip ^U {ia, xe'^ 4 ) 



W(a,-x) 



k(a 



2 e i™9 



(4.10) 



a) u{ 



ia, xe 



These relations are convenient for numerical computations because 
for x > and x < we can use the same [/—function. 

4-1.1 The function p(a) We give more details on the function p 
defined in (4.9). For large values of a it is convenient to use the 
representation 



p(a) = \tt - \a + \a\xia 2 + p*(a), 



(4.11) 



where p*(a) = 0(1/ a) as a — > oo. To give more details we give an 
asymptotic expansion. We have Binet's formula (see [9], p. 55, for an 
integrated version) 



In r(z + J) = z ln(i + z)-^-z + i ln(27r) 



(4.12) 



lo 



where 



(4.13) 



with Cfc in terms of Bernoulli polynomials: 



t\ < 2vr, 



fc fc + 3 
2^+2" 



, fc = 0,l,2,. 



This gives the asymptotic expansion 
as ±a — > oo, where 



<4 = (-1)" 



C2fc 



(2fc + l)(2fe + 2)' 



fc = 0,1,2,... 



(4.14) 



(4.15) 



(4.16) 
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The first few coefficients are 

j _1_ ,7 = 13. j 37 j 29 j 1129 o i ? \ 

"0 12 ) "1 720' " 2 20160' 3 26880' 4 1520640' 



4-2 Integral representations for a < 

For W{— a,±x) we consider (2.1) for U ia, xe"- 7 ™/ 4 ) (see (4.7) and 
(4.10)), that is, 

where C is a vertical line on which 3f?s > 0. On C we have — \-k < 

phs < ^7r, and the many- valued function s ia ~2 assumes its principal 
value. The transformations 

x = 2t^/a, s = yfaw (4.19) 

give 

U (-ia,xe-^) = f 

V / iy/2n Jc Vw 

where 

(j)( w ) = l w 2 - 2te~ ni/4 w + i\nw. (4.21) 
The saddle points follow from solving 



W - 2te~' Ki / A w + i 



#W = ^ ™ w - ' - = 0, (4.22) 
giving saddle points 

w±=u± + iv± = e~ ni/A (t ± Vt 2 + l) , v± = -u±. (4.23) 
The relevant saddle point is w+. We have 

( j ) { w+ ) = \ 1T + l (t 2 + 2i-\), (4.24) 

where £ is given in (2.12). The path of steepest descent through w + 
is for \6\ < \tx defined by 

±r 2 sin2# - 2trsm(6 - W) + lnr = = t 2 + 2£ - \, (4.25) 
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where w = re ld . In rectangular coordinates w = u + iv this equation 
reads 



uv 



+ V2t(u - v) + \ ln(u 2 + v 2 ) = i 2 + 2£ - ±. (4.26) 



We can solve equation (4.25) for sin(# — \tt) (it is a quadratic equation 
for this quantity), giving 9 as function of r. This makes it possible 
to integrate (4.20) with respect to r, but this introduces singularities 
in the integral where r attains its minimal value, although the path 
itself is smooth. 

Integrating with respect to 8 or v is a better option. We can nu- 
merically determine the path in an algorithm, but this is not a very 
efficient method. Instead, we replace the steepest descent path de- 
fined in (4.26) by a path u(v) such that 

1. u(v) is smooth for all v £ IR; 

2. u(v) passes through the saddle point: u(v + ) = u + ; 

3. du/dv at v + has the same value as du/dv for the steepest descent 
contour at w + ; 

4. the path runs into the valleys of e^ w ^ at ±zoo. 

^From (4.26) we can show that du/dv = at the saddle point w+. 
Hence, a simple path C that fulfills the four conditions is the vertical 
line u = u + . Introducing q = v — v+, using w = w+ + iq and (4.24) 
we obtain for (4.20) the representation 



U (y—ia, 



xe 



/oo 
e- a ^g(q)dq. 
-oo 



(4.27) 



The function ip(q) is given by 



</>(<?) = <j>{w+) - 4>(w), g(q) = . (4.28) 

VI + zq/w + 



For small values of q we have 



m = 1J >^ + ^+o{q% 



(4.29) 



n + = ^(t + v^TT). 
We conclude that W(—a,±x) (see (4.7) - (4.11)) are given by 



V 71 " F+l 



/oo 
e- a ^g(q) dq 
-oo 



(4.30) 
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W(-a,-x) 



i 

at 



^tt \w+\k(-a) 



f c 

' L 



-aip(q) 



g(q) dq 



where 



(4.31) 
(4.32) 



For the derivatives we find, starting with (4.18), 



V A U' (-m,xe"™/ 4 ) 



e \-na-\ia a \ + \ai 



f 

J — i 



-aip{q) 



h(q) dq, 



where 



h{ q )=(y¥Tl-e-^' A q)g{ q ). 
It follows from (4.10) that W'(—a,±x) are given by 



i e 



>X 



/oo 
e~ a ^i)h(q)dq 
-oo 



(4.33) 



(4.34) 



(4.35) 



3 

a* 



/oo 
e -^(9)h(q) dq 
-OO 



(4.36) 



W r '(-a,-x) = ^ /< 

V 7T |w+|fc(-a) 

For large values of a and/or t the oscillatory behaviour of W(— a, ±x) 

and W'(—a, ±x) is mainly described by the exponential factor e 2ia ^ 
contained in e %x . The other elements of these formulas are slowly 
varying. 

Asymptotic expansions follow from [10]. Formula (2.29) of that 
paper gives, with fj, = \/2ae~ m ^, 



a*(t 2 + 1)3 



(4.37) 



where 4> s are polynomials given in (2.11) and r in (2.32) of [10]. 
Formula (2.33) of [10] gives 



W'(-a,x) ~ ^Jk(-a)a^(t 2 + 1)1 3? 



° H)W) 



?, e 



''X 



E 

s=0 



(2a)* 



, (4.38) 



where V>s are polynomials given in (2.16) of [10]. 
For W(— a, —x) and its derivative we have 



W(— a, —x) ~ — 



a*(t 2 + 1)3 



s=0 



(2a) 



s 



(4.39) 
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and 

W'(-a,-x) ~ - 



k(-a)ai(t 2 + l)4 9 



oo 
s=0 



(~0 s ^(r) 
(2a)- 



, (4.40) 



The asymptotic expansions in (4.37) - (4.40) hold when a — > oo, 
uniformly with respect to t > —t , but also for t — > oo, uniformly 
with respect to a > ao, where ao and to are fixed positive numbers. 



4-3 Integral representations for a > 

Because of the turning points we consider three cases. We write x = 
2ty/a. We use the [/—function in (4.7), and write (2.1) in the form 

U (ia, xe'^) = 4^ I e -^ i/4s+ ^s- ia ^=, (4.41) 
with conditions as in (4.18). The transformation s = yfaw gives 

U Ua,xe~^) = <L±^t^ f (4 . 42) 



where 

(f){w) = l w 2 - 2te~ wi/4 w -ihxw. (4.43) 



4-3.1 The case t > 1 The saddle points are now 

w± = u±+ iv± = e~ ni/4 (t ± Vt 2 - l) , v± = -u±. (4.44) 

The relevant saddle point is w+ , and for numerical integration a con- 
venient choice of C is the vertical line through w+ . 
Using 

-\ix 2 + a(j)(w + ) = 2ia£ + \ia - \na, (4.45) 

where £ is given in (3.38), and writing in (4.42) w = w + + iq, we 
obtain the analogue of (4.27) 

U (ia, xe~™l A ) = J e 2ia? J e~ a ^ q) g{q) dq, (4.46) 

where 

iP(q) = cP(w + )-cP(w), g{q)= _ L= . (4.47) 

y/1 + iq/w+ 
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It follows that 



W(a, x) = Vpp^ K [F(o> x)] (4 . 48) 



and 



where 



W(a, -x) = V^°l° 4 Q [F(a, x)] . (4.49) 



where 



(4.51) 



1 /*oo 

F(a,x) =e i [p*W+4-+ 2i <] / e - a ^g(q)dq (4.50) 

J — oo 

For the derivative we find, as in (4.33), 

/oo 
e - av,(,?) %)dg, 
-oo 

where 

%) = (e-^A/P^l + iq) g(q). (4.52) 

It follows that 

W'{a, x) = - ^y°, 4 5ft [G(a, x)] (4.53) 

y/Tr\W + \ 

and 

3 

W(a, -x) = J' 9 [G(a, x)] . (4.54) 



G(a,x) = e *[p*W+4-+ 2i <] / e - a ^ q) h(q)dq (4.55) 

J— oo 



Asymptotic expansions follow from (2.9) of [10]. By changing \i 
fie-™/* in that formula we obtain 



a* (i 2 - 1)4 



i[p*(a)+i7r+2<] y j^sCO 



, (4.56) 



where </> s (t) are the same polynomials as in (4.37), and 

1 / t 



2 VV^T 



1 . (4.57) 
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Formula (2.18) of [10] gives 

W'(a,x) ~ -^/k{a)a^(t 2 - 1)2 



E 

s=0 



(2a)- 



(4.58) 



where ip s ( T ) are the same as in (4.38). 

The asymptotic expansions in (4.56) and (4.58) hold when a — ► oo, 
uniformly with respect to i > 1 + to, but also for t — > oo, uniformly 
with respect to a > ao, where ao and to are positive numbers. 



4.3.2 The case -1 < t < 1 We use (4.41), (4.42) and (4.43) with 
saddle points 

«, ± = e -W4 ^ ± i^l-t 2 ) = e-Wfct* t = cos (4.59) 

which are located on the unit circle. We have 

(j>{w ± )=±2r ] -\^ + i{\+t 2 ), (4.60) 

where rj = \{6— sin cos 9) is also used in §3.1 and defined in (3.9). We 
see that the imaginary parts of <f>{w±) are equal. As a consequence, 
the steepest descent path may go (and in fact in the present case does 
go) through both saddle points. 

In Figure 3 we have shown the paths for three values of t. The 
contours run from — ioo to W-, then along the arc to w+ (in the 
direction of the arrows), and from w+ to +ioo. Through each saddle 
point the local contours of steepest ascent and steepest ascent are 
shown. The complete contours include steepest descent parts and 
steepest ascent parts. 

^From (4.60) we see that tu+ is dominant for < t < 1 (rj is 
positive for these values of t). Another point of interest is that the 
oscillatory factor e~± lx2 = e _m * 2 in (4.42) is nullified when we put 
4>(w-) or (f>(w + ) in front of the integral, because ^s<fi(w±) = \ + 
t 2 . This explains that the function W(a,x) does not oscillate if t € 

[-1,1]- 

When we write w = u + iv and integrate in (4.42) with respect to 
v we obtain 



U (ia, xe 

PCX) 

/ e " 

JV- 



) 



1 

a*" 



■\ai e a\\i- 



2vr 



■ai>i(v) 



h(v) dv + e- 4ar > [ V e- a ^f 2 (v)dv 

J —oo 



(4.61) 
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where 



tpi(v) = - <p{w), ip 2 (v) = 4>{w-) - 4>(w), 



(4.62) 



and 



fi(v) 



,du\ 
> f 

dv 



.du 2 
dv 



(4.63) 



We may assume different relations between u and v in both integrals; 
this explains u\ and u 2 , which are functions of v. 




Figure 3. Steepest descent contour for the integral in (4.61) 
for several t— values. The contours run from — zoo to u>_, then 
along the arc to u> + (in the direction of the arrow) , and from w + 
to +ioo. Through each saddle point the local contours of steep- 
est ascent and steepest ascent are shown. The complete con- 
tours include steepest descent parts and steepest ascent parts. 
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It follows from (4.10) that 



W(a,x) 



W(a,-x) 



ate 



7T 



(4.64) 



a 4 e 



s T| K(a,x) 



where K(a, x) denotes the sum of the integrals between the square 
brackets in (4.61). 

For the derivatives we find, starting with (4.42), 

_I,-,.2 31, 



U' (ia,xe- ni ^) = 
and using (4.10), 



2vr 



■ jf e a ^ w) (te 



-ni/4: _ 



W 



dw 



(4.65) 



W^x)=^ a } e2a \ 



jWW-fr} K d {a,x) 



W'(a, -x) 



ah 2ar] 



(4.66) 



where Kd(a,x) denotes the sum of the integrals between the square 
brackets in (4.61) with fi(v) replaced with gi(v) = (te~± nt — w)fi(v). 

The asymptotic expansions follows from (2.29) and (2.33) of [10]. 
We change t — > —it and \i — > /Lie^ 7 ™ in these formulae and obtain 



m, xe 



-7ri/4 N 



E 



U' (ia, xe"^/ 4 ) 



a -^ a +i e 5 ia+ l^-i™ +2OT '(l - t 2 )i (4.67) 



E 



X j> z 

s 



V2 



o (2«) s ' 



where r\ is defined in (3.9) and 4> s {t*) are the same as in (3.35). 
It follows from (4.10) that 



W(a,x) 



(24 



-e 2 ^ 



W"(a,s) ~ -^/fc(a)a*(l - t 2 )\e 2ar > 



(4.68) 



o W 
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The asymptotic expansions in (4.68) hold when a — > oo, uniformly 
with respect to t € [— 1 + S, 1 — S], where 5 is a fixed positive number. 

4-3.3 Unstable representations For large values of ar\ the representa- 
tions for W(a, —x) and W'(a, —x) in (4.64) and (4.66) are unstable. 
To see this, observe that (4.68) can be used for t £ [— 1 + 5, 1 — 5]. 
The dominant behaviour comes from \fk{a) e 2av . Since (see (4.8)) 
k(a) ~ ^e~ na , the dominant behaviour comes from e ax , where x = 
arcsin t—ty/l — t 2 , an odd function that is positive on (0, 1] . This dom- 
inant behaviour does not appear in the representations for W(a, —x) 
and W'(a, —x) in (4.64) and (4.66). There we see the dominant parts 
e 2ari I y / k(a); in 77 we use positive t when x is positive. It follows that 
the imaginary parts in the right-hand sides of (4.64) and (4.66) have 
to be very small when an is large. In fact, the first integral in (4.61) 
should be of order e~ 4av in that case, which is not apparent from this 
representation. 

A possible solution to this problem is using the representations for 
W(a,x) and W'(a,x) in (4.64) and (4.66) for t £ [-1,0]. However, 
when t i — 1 the phase of w + becomes 3tt/4 and that of becomes 
— 5ir/4, which is outside the standard interval (— it, tt] of the phase 
of w in (4.42); that is, w is outside the standard Riemann sheet. In 
Figure 3 the path for the case t = —0.99 is shown. The technical 
details will be worked out when writing the numerical algorithms. 

5 Concluding remarks 

In a future paper we will discuss the numerical aspects and describe 
computer algorithms based on the integral representations given in 
this paper. Several quantities have to be calculated with great care. 
For example, straightforward use of ip(9) defined in (2.10) when 9 
is small, that is, at the saddle point, will give cancellation of lead- 
ing digits. Also, to represent the functions for a large range of the 
parameters scaling is needed. 

When implementing the representations we will decide if the steep- 
est descent paths will be used or approximations of these paths, as 
we suggested for the W— function in § 4.2. For example, integrating 
in (2.3) along the vertical line through the saddle point wq gives a 
simpler representation than (2.9). However, the integral along the 
vertical line has a non-real phase function. Another approximation of 
a steepest descent contour is given in (3.13). We will investigate effi- 
ciency aspects in combination with programming aspects in deciding 
which representation in these examples should be used. 
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This also holds for the quite complicated steepest descent paths 
in § 4.3.2. We have not indicated in (4.61) the relation between v and 
u on the different parts of the path. This will be done during the 
implementation of the algorithms. 
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